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Abstract 



We propose a class of semi-Lagrangian methods of high approximation order in 
space and time, based on spectral element space discretizations and exponential 
integrators of Runge-Kutta type. The methods were presented in [7] for simpler 
convection-diffusion equations. We discuss the extension of these methods to 
the Navier-Stokes equations, and their implementation using projections. Semi- 
Lagrangian methods up to order three are implemented and tested on various 
examples. The good performance of the methods for convection-dominated 
problems is demonstrated with numerical experiments. 
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1. Introduction 

Consider the incompressible Navier-Stokes equations 



Here u = u(x, t) is the velocity field defined on the cylinder Q. x [0, T] {O. c R"^ 
for d = 2, 3), subject to the incompressibility constraint ([2]), while p = p(x, t) is 
the pressure and plays the role of a Lagrange multiplier, and v is the kinematic 
viscosity of the fluid. We consider no slip, or periodic boundary conditions when 
the domains allow it. 
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Uf -I- u Vu = vV^u - Vp, 
V • u = 0. 
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For no slip boundary conditions we will mostly consider the case 



Vihn = 0. 



(3) 



The variables (u,p) are sometimes called primitive variables and the accurate 
approximation of both these variables is desirable in numerical simulations. 

A typical approach for solving numerically convection-diffusion problems is 
to treat convection and diffusion separately, the diffusion with an implicit approach 
and the convection with an explicit integrator, see for example El [11 EZl. We 
will refer to these methods as implicit-explicit methods (IMEX). The advantage 
of this approach is that most of the spatial discretizations of the diffusion operator 
give rise to finite dimensional counterparts which are symmetric and positive 
definite, so the implicit integration of the diffusion requires only the solution of 
symmetric positive definite linear algebraic systems. 

In this paper we propose high order semi-Lagrangian discretization methods 
in time, to be used in combination with high order spatial discretizations of 
the Navier-Stokes equations, as for example spectral element methods. High 
order methods are particularly interesting when highly accurate numerical 
approximations of a given flow are required. An interesting field of application 
is the direct numerical simulation of turbulence phenomena, as pointed out for 
example in [|36l . Another relevant situation is in connection with discontinuous- 
Galerkin methods, as an alternative to the use of explicit Runge-Kutta schemes. 
In this context, the purpose is to alleviate severe time-step restrictions imposed 
by CFL conditions. See for example (26l for the use of IMEX time-stepping 
schemes combined with discontinuous-Galerkin space discretizations, and [l32ll 
for semi-Lagrangian discontinuous Galerkin methods. 

The integration methods we propose in this paper are implicit-explicit 
exponential integrators of Runge-Kutta type. They combine the use of 
a diagonally implicit Runge-Kutta method (DIRK) for the diffusion with a 
commutator-free exponential integrator (CF) for the convection, and were denoted 
DIRK-CF in [[6l|71. In addition to being implicit-explicit, because of the presence 
of exponentials of the linearized convection, the methods are amenable for semi- 
Lagrangian implementations providing improved performance in convection- 
dominated problems. In the present paper we extend the approach of [7| to the 
incompressible Navier-Stokes equations and we assess the performance of the 
obtained semi-Lagrangian schemes. 

To clarify the connection between semi-Lagrangian methods and exponential 
integrators, let us consider the simple linear convection-diffusion model problem 



Du 
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where ^ := ^ + V ■ V is the total derivative of u, with V : R'' ^ R"^ 
the convecting vector field. Semi-Lagrangian methods for convection-diffusion 
problems are usually formulated by applying a time-discretization directly to the 
time-dependent PDE problem before any space discretization is performed. A 
simple example of a semi-Lagrangian scheme is 

Un+l — Un(X(t„)) 2 
= yV Un+l, 

h 

where Kit) is the solution of 

Kit) = \iXit)), Xitn^i) = X, X 6 c R^. 

This amounts to discretize the total derivative with forward Euler and the diffusion 
with backward Euler. Finally the space discretization is introduced. The practical 
realization of this method requires: 

• introducing a space discretization, where U„ is a vector representing the 
numerical solution at time t„ and on all nodes of the discretization grid F (or 
Un belonging to a suitable finite element space); 

• a discretization of the Laplacian; 

• an operator J(/„ interpolating t/„ on the domain Q; 

• a suitable integration method to solve the equation X = Y(X), and compute 
the characteristic paths; we denote by 0/,(x) its numerical flow at time h and 
with initial value x e F. 

So the fully discrete method can be expressed in the form 

- J[/„(0,(F)) = h vAUn^i. (4) 

Our approach is to formulate the semi-Lagrangian methods as time-integrators 
of appropriate semi-discrete PDE problems. For this purpose, we assume that 
there exists a discrete convection operator C(V) (in matrix form), such that 

Jf/(0„(F)) = exp(/zC(V)) U, (5) 

where exp denotes the matrix-exponential, and V is a discrete version of V and 
is typically known only on the nodes of the discretization grid F. Based on this 
assumption the semi-Lagrangian method Q can be also written in the form 

Un^i - exp(/zC(y)) Un = hvAU„^,. 



3 



This corresponds to the time-integration of a system of ordinary differential 
equations (ODEs) of the type 



U{t) + C(V)U(t)=AU(t), (6) 

with an exponential integrator. As a result the semi-Lagrangian methods of [3 
can be conveniently formulated as exponential integrators for a semi-discrete 
convection-diffusion problems of the type ([6]) and (nonlinear variants of ([6])). 

In the case of pure convection problems a well known estimate, due to Falcone 
and Ferretti [[Ml [151, gives a bound for the local error t(x„ f„+i) of the type 

\T(Xi,tn+l)\<K\h'' + 

where is a generic grid-point and tn+i is time. In this estimate, the term h'' is 
the error due to the numerical approximation of the characteristics paths, the term 
arises from the accumulation of the interpolation error, and ^ is a constant 
independent of h and Ax. This estimate of the error suggests that the spatial error 
is affected positively by the use of large time-steps h, moreover, when high order 
interpolation is used (like in the case of high order spectral element methods), the 
integration of the characteristic paths should also be done at high accuracy and in 
particular an optimal h could be chosen so that 

h 

This motivates the interest in designing semi-Lagrangian methods achieving high 
order in time, when the adopted space-discretization is of high order. In [7J and 
[[9l we considered high order space discretizations by spectral element methods 
for convection-diffusion problems and provided high order semi-Lagrangian time- 
discretizations for nonlinear convection-diffusion problems. We also showed 
numerically that the proposed integrators do overcome nominal CFL stability 
restrictions. 

So far the case of linear and nonlinear convection-diffusion equations have 
been considered. Semi-discretizations of the Navier-Stokes equations, giving rise 
to index 2 differential-algebraic systems, have been approached successfully by 



Notice however that this underlying operator C described above will never be used explicitly 
in the implementation of the semi-Lagrangian methods, and the exponentials of C will be obtained 
computing characteristics. With abuse of notation, we will also denote by C the semi-discrete 
convection operator based on spectral element space-discretizations, which will be used for 
comparison methods based on classical (Eulerian) IMEX time-integration. 
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BDF-like multi-step methods proposed in [9]. The connections of these methods 
to the methods proposed in OOll . [[36l and [1211 . have also been explained. 

Here we address the case of semi-Lagrangian methods based on one- step 
formulae, and more precisely Runge-Kutta type formulae. Given a time-stepping 
technique, a standard approach to adapt the method to the incompressible Navier- 
Stokes equations is by means of projections. The primary example of this 
technique, and most famous projection method for the incompressible Navier- 
Stokes equations is the Chorin's projection method, proposed by Chorin [|TOHTT]| 
and Temam [35]. As for the usual formulation of semi-Lagrangian methods, 
this and similar projection methods are also typically formulated by applying a 
time-discretization directly to the time-dependent PDE problem before any space 
discretization is performed. The methods appear as splitting methods, where 
solutions of unconstrained problems and Poisson equations for the pressure are 
handled separately and composed. 

The study of the temporal order of the Chorin's projection method was 
considered in [l34ll and [31] and it revealed order 1 in time for the velocity and only 
^ for the pressure. Such loss of order in time for the pressure has been analysed 
for projection methods for Stokes and Navier-Stokes equations, |[23]| . and there 
are various ways to add correction terms and restore the full order. 

We instead choose a different strategy. We first semi-discretize in space, 
taking care of boundary conditions, and then apply the exponential integrators and 
perform projections at the space-discrete level. To define the appropriate semi- 
discretization of the convection, we make use of the assumption (|5]). We start with 
a time-dependent spectral element space-representation of the numerical solution 
taking care of boundary conditions. The space-discrete numerical approximations 
of the Navier-Stokes velocity and pressure solve a system of index 2 algebraic 
equations with linear constraints. These semi-discrete equations are described 
by the usual semi-discrete divergence operator, a corresponding Lagrangian 
multiplier representing the discrete pressure, and the usual discretization of the 
linear diffusion operator. The discrete nonlinear convection C{U) is assumed to 
have the property that for any two approximations of the numerical solution U 
and Un, ^ holds with V = U and U = Un. 

We then eliminate the Lagrangian multiplier from the discrete equations 
and apply the semi-Lagrangian exponential integrators to the resulting system 
of unconstrained ordinary differential equations, ensuring incompressibility in 



two different ways described in section 3 The numerical approximation of the 



pressure is obtained in a post-processing step. 

In section 2 we discuss various projection techniques for Runge-Kutta splitting 
methods (IMEX and semi-Lagrangian methods) applied to the Navier-Stokes 
equations, to motivate our approach. In [section 3] we present our schemes and 
show how to obtain high order implicit-explicit and semi-Lagrangian methods 
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for the incompressible Navier-Stokes equations, section 4 is devoted to the 
description of the implementation details. In this section we describe how we 
use techniques from [|17ll for the efficient solution of the linear algebraic systems. 
We obtain an overall strategy which resembles conventional projection schemes 
as described in [23 J, where, at each step in time, one only needs to solve a 
sequence of decoupled elliptic equations for the velocity and the pressure. In 



section 5 we report the numerical experiments. We provide numerical verification 
of the temporal order of the methods, and we demonstrate the clear benefits 
of the proposed semi-Lagrangian methods in the case of convection-dominated 
problems. 

2. Projection methods for the incompressible Navier-Stokes equations 

Since the pressure in the Navier-Stokes equations is not a dynamic variable, 
but plays the role of a Lagrangian multiplier, a usual approach to time integration 
of the equations is by diff"erentiating the constraints with respect to time until a 
diff"erential equation for p is obtained. 

Another option is that of eliminating p from the equation ([T]). This can be done 
by taking the divergence of the equation ([T]) (equivalent to differentiating ([2]) once 
with respect to time). This way we obtain a Poisson equation for the pressure 

= -V-(u- Vu), (7) 

with boundary conditions to be deduced from 

u^n + (u • Vu - vV^u)! = - {Vp%^ . (8) 

In particular, when u is space-periodic, the pressure p is fully defined in terms 
of the velocity field u and the periodicity condition. In the case of no slip 
boundary conditions, taking the inner product of ([8]) with the unit normal n gives 
the Neumann boundary conditions: 

^ = yV^u ■ n, on 
an 

These can be used to solve the Poisson equation for p and fully determine the 
pressure. In both cases we can write p = i/r(u) (see [fT9l Ch. II]). We can then 
eliminate the pressure from ([T]) and obtain 

u, = yV^u - u • V u - Vi/r(u). (9) 

At this point, (|9]) contains only dynamic variables and time integration 
methods can be applied in a usual way. We want to discuss the use of a time 
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integration method on this equation and the preservation of incompressibility 
under time discretization. Suppose we apply a classical Runge-Kutta method to 
this equation, without assuming any discretization in space. Provided the initial 
value is divergence-free, since the update of a Runge-Kutta method is always of 
the form 

.s 

i=l 

where Z?,- for / = 1 , . . . , 5 are the weights of the Runge-Kutta method, and 

F; = (vV^u - u • Vu - V(A(u))| , U; « u, i=l,...,s, 

and since V • F, = 0, then u„+i is also divergence-free. The order of the numerical 
approximation of the velocity in time is the order of the Runge-Kutta method (see 
[[^ for a similar conclusion). An approximation for the pressure of the same 
order can be obtained by p„+i = (/r(u„+i), i.e. solving (|7]) with u„+i inserted at the 
right hand side and with appropriate boundary conditions. This approach gives 
high order divergence-free approximations of the velocity, but it lacks flexibility 
because convection and diffusion are treated with the same Runge-Kutta method, 
either an explicit or an implicit one. On one hand the use of explicit Runge-Kutta 
methods usually leads to severe step-size restrictions due to the presence of the 
Laplace operator, but on the other hand the use of implicit methods involves the 
solution of a large nonlinear system of equations at each step. A more popular 
approach is therefore to treat convection explicitly and diff"usion implicitly by a 
splitting and composition method or, alternatively, by a so called IMEX method. 

If we use a splitting and composition method to solve (|9]) numerically, then we 
need to make sure that the individual flows which are composed are divergence- 
free in order for the final approximation to have the same property. Suppose we 
can write if/(u) = ij/i{u) + 1/^2 (u) such that we can split the right hand side of (|9]) 
into two divergence-free parts 

/i(u) = yV^u - V<Ai(u), /2(u) = -u ■ Vu - V^2(u), (10) 

if we then solve individually each of the two split parts (either exactly or by 
applying a method guaranteeing the incompressibility of the respective numerical 
approximations), then also the composition of the two numerical flows will 
be incompressible. The two terms (/'i(u) and i^2(u) arise as solutions of two 
corresponding Poisson problems imposing appropriate boundary conditions. A 
difliculty when performing such a splitting of the PDE is that these boundary 
conditions should be consistent with those imposed for i/r(u) in (|7]), and such that 
the two individual problems Uf = /(u) with i = 1, 2 in ( fTO] ) admit solutions. In the 
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next sections we deal with this issue by first discretizing in space ([T])-([2]) with a 
spectral element discretization, then eliminating the discrete pressure Lagrangian 
multiplier and finally applying the integration methods to the discrete analog of 

u, = yV^u-u-Vu-V^i(u)-Vi/r2(u), (11) 

see section [3n 

The situation for semi-Lagrangian and IMEX methods is similar to what 
happens with splitting methods. We clarify this analogy in two simple examples. 

Semi-Lagrangian methods can be seen in some cases as splitting methods. If 
we consider the reformulation of ([9]) using the total derivative we obtain 

^ = vV2u-V.A(u), (12) 

and a semi-Lagrangian method is based on the direct approximation of the total 
derivative along characteristics. So for example using a simple forward Euler 
scheme we get 

U„+i =Un + hvV^Un - V(A(u„), u„ := Un(X(t„)) (13) 

where X{t) are the characteristic paths, and the convecting vector field is V := u„. 
Now u„ = Un(X{t„)) = u{h) is the solution at time t = h of the pure convection 
problem 

Du 

— = 0, u(0) = u„, 

and u„ is not necessarily divergence-free even if V ■ u„ = 0. The given semi- 
Lagrangian method is evidently a splitting and composition method where we 
compose two numerical flows which are not divergence-free, then we can not 



expect u„+i to be divergence-free. In section 3.3 we consider splitting the discrete 



counterpart of the Navier-Stokes equations in the form (11), and we apply different 



methods to the discrete versions of pO) ), with this strategy we obtain an adaptation 
of the DIRK-CF methods of [7J to the Navier-Stokes equations. 

The situation is similar for IMEX methods. The simplest IMEX method is 
obtained by combining the forward Euler and backward Euler method and it is 
also a splitting method. Advancing the numerical solution of (|9]) by first applying 
a forward Euler step to the convection, and then a backward Euler step to the 
diffusion and the term ij/(u), we get 

u„+i = u„ - /?u„ ■ Vu„ -h /?vV^u„+i - /zVi/f(u„+i), 

and with p„+i = (/r(u„+i) the approximation is of order 1 in time, both for the 
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pressure and for the velocity, but V ■ u„+i 0. If we instead use the method on 



( [TT] ) then we obtain 

U„+i = U„ - /ZU„ • Vu„ - /zVl/r2(u„) + /zvV^U^+i - /lV(/ri(u„+i). (14) 

In this case u„+i is a sum of three divergence-free terms and is therefore 
divergence-free. This method has order 1 for the velocity, and with = i/r(u„+i) 
we also get an order 1 approximation for the pressure. We will pursue this strategy 
for IMEX methods in the next sections. 

We want to emphasize the difference between this approach and the Chorin's 
projection method which consists of two steps: an update step 

u„+i = u„ - /zu„ • Vu„ -h /zvV^u„+i, (15) 

and the projection step 

u„+i = u„+i -h V A;?„+i = -V • u„+i, (16) 

guaranteeing incompressibility of u„+i (and we omitted boundary conditions for 
for simplicity). In the Chorin's projection method one Poisson problem for 



the pressure is solved in ( [16] ). In this case the pressure update does not 
satisfy Q with u„+i inserted at the right hand side. As shown in [l34l[3T]| . pn+\ is 
an approximation of order \ for the pressure. 

We finally remark that given a differential equation whose solution is evolving 
on a manifold, a simple approach to obtain a numerical approximation of high 
order evolving on the manifold is by applying any one- step integrator to the 
problem, and then performing an appropaiate projection on the manifold (see 
AM IV.4], see also Sect.5.3.3] and [|25l Sect.VII.2]). For the Navier- 
Stokes equations the solution is constrained to evolve on the linear subspace of 
divergence-free vector fields. So any integration method applied to (|9]) giving an 
approximation of high order for the velocity, can be corrected at the end of each 
step by a suitable projection step to enforce the incompressibility constraint. A 
final post-processing step, projecting the right hand side of the equation evaluated 
at the obtained approximation of the velocity, leads finally to an approximation of 



the same order for the pressure. This approach will be pursued in section 3.4 



3. High order implicit-explicit and semi-Lagrangian methods of Runge- 
Kutta type for the incompressible Navier-Stokes equations 

In this section we present the details of the high order integration schemes 
proposed in this paper. 
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3.1. Semi-discretization and elimination of the discrete pressure variable 

As outlined in the introduction and in the previous section we consider the 
spatial discretization of ([T]) and proceed eliminating the Lagrangian multiplier 
from the semi-discrete equations. 

After spatial discretizations of type spectral-Galerkin or spectral element 
methods, we obtain a system of differential-algebraic equations of the type: 

By = Ay + C(y)y-Dh, 
Dy = 0. 

Here the matrix A represents the discrete Laplacian, C(y) is the discrete convection 
operator; D and are respectively the discrete divergence and gradient 
operators. Meanwhile vectors y and z represent numerical approximations of the 
velocity u and the pressure p respectively. In the numerical experiments A and D 
are obtained by a spectral element method (SEM) based on the standard Galerkin 
weak formulation. 

For C(y) we have two options. When the overall space-time discretization 
is of Eulerian type, and we integrate in time with IMEX Runge-Kutta methods 



like ( [T4| ) (see also Appendix A.l I, then C{y) is assumed to be the semi-discrete 
convection operator based on the standard spectral element space-discretization. 
Otherwise, when the overall space-time discretization is of semi-Lagrangian type. 



and we integrate in time with DIRK-CF exponential integrators like ( 13 1 (see also 



[Appendix A.l ), then C(y) is the discrete convection operator fulfilling ([5]), and 



Jj/ is the natural, built-in interpolation operator of the SEM method. 

In the numerical experiments, the approximation is done in - Vn~2 
compatible velocity-pressure discrete spaces. That is, in each element we 
approximate the velocity by a N-degree Lagrange polynomial based on Gauss- 
Lobatto-Legendre (GLL) nodes in each spatial coordinate, and the pressure by 
(A'^ - 2)-degree Lagrange polynomial based on Gauss-Legendre (GL) nodes. The 
discrete spaces are spanned by tensor product polynomial basis functions. A 
consequence of this choice of the spatial discretization is that the grid for the 
pressure does not include boundary nodes (there are no boundary conditions 
for the pressure). We note that this is not the only viable choice of spatial 
discretization for our time-integration schemes, but, for the sake of clarity, we 
find it preferable to fix a concrete spatial discretization. 

We can eliminate the Lagrangian multiplier z from the system ( [T7] ) by 
multiplying on both sides by and using the constraint Dy = 0. This leads 
to the following system of ordinary differential equations (ODEs) 



By = Ay + C(y)y-BHB-\A + C(y))y, (18) 
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where H .= B ^D^(DB ^D^) ^D. We introduce the projection 

Ii:=I-H 



allowing to write the ODE in the short form 

y = IiB'^Ay + IlB'^C(y)y. (19) 

3.2. IMEX methods applied to the projected semi-discretized Navier-Stokes 
equations 

To begin with we apply the IMEX Runge-Kutta method to the system and 
we obtain an high order implicit-explicit Eulerian approach for the incompressible 
Navier-Stokes equations. 

We assume that the considered IMEX methods are stiffly accurate, i.e. the last 
stage coincides with the update of the solution. This assumption is not strictly 
necessary, but since in this case the update j^+i coincides with the last stage value 
of the method 7^, the format of the methods and the formulae are simpler. We 
obtain 

tov i= \ : s do 

Yi = yn + hliB ^ 2 (aijA + aijC{Yj)) Yj + haaUB-^AYi 

7=1 ^ ^ 

end for 

yn+i = Y,. 



See Appendix A.l for the definition of IMEX methods for convection- 
diff"usion problems. The method requires the solution of one linear system per 
stage with the matrix (I-haijIlB'^A). This matrix is never formed explicitly in the 
implementation. Instead the equation can be rearranged to obtain a method for the 
diff'erential- algebraic equation ( [T7| ) and, denoting by Z/ := {DB~^D^)~^DB~^A 7/, 
/ = I, . . . , s, we obtain 



(-1 



BYt - ha^AYi + hai^iD^Zi = By„ + hJ^ (cHjiAYj + D^Zj) + ai,j{CiYj)Yj - D^Zd) 



7=1 



DYi = 0. 



We notice in particular that = 7, satisfies the discrete incompressibility 
constraint Dyn+i = 0. 
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Finally, to recover the correct approximation of the pressure at time t„+i we 



perform a post-processing step. We consider the right hand side of ( [17] ) and 
evaluate it in y„+i obtaining an approximation to Dy(tn+i). Since Dy„^i = and 
Dy{tn+\) = the correct approximation of the pressure is given by Zn+i such that 

B-'Ayn^i + B''C(yn+,)yn-,i - 5"'Z)^z„+i = 0, 

and amounts to solving a linear system for Zn+i, obtained by multiplying by D: 

DB-'D^Zn+i = DB-\A + C{yn^,))y„^,. (20) 

3.3. DIRK-CF methods applied to the projected semi-discretized Navier-Stokes 
equations 

In this section we present our first semi-Lagrangian approach to the 
incompressible Navier-Stokes equations based on the exponential integrators 
methods of [lZJ. 



We apply the DIRK-CF methods to the equations ([19]). See Appendix A.l for 
the definition of DIRK-CF methods. We get 

for z = 1 : 5^ -I- 1 do 

Yi = <piy„ + h Z aijifiif^'UB-'AYj + hauUB-'AYi 

7=1 

Yj:= Z.<i'.forr=l,...,7 
(Pi = expihUB-^C(Yf)) ■ ■ ■ exp(hUB-^C(Yl)) 
end for 

yn+i = Ys+l, 

where a^+ij = bj, aj^j ^ = fi'^, j,k = I, . . . , s, and a^+i = 0. Also in this 
case the methods can be reformulated as methods for the differential-algebraic 



equation ( [T7[ ) as done in the previous section. We will not repeat this here and 
observe instead that also in this case the stages satisfy DY, = 0. The exponential 
exp(hnB~^C(w)) ■ g is the solution of the semidiscretized equation 

Bv = Ciw)v + D'^z, (21) 
Dv = 0, (22) 

on the time interval [0, h]. This differential equation is the discrete counterpart of 
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a system of linearized Euler equations which we can write in the form 

y^ + \.Vy = Vp, y{xuQ)=gi, in [Q,h]x^, 
V ■ 7 = 0. 

To approximate each of the exponentials exp(/zn5~^C(w)) • g we envisage two 
options: 

• Use a projection method of high order for 7: 

exp(/in5-^C(w)) • g = UexpihB~^C(w)) -g + UE. (24) 

• Consider the vorticity formulation of 0}f + \- Vo) + t{o)) = and o) = 
V X 7, and approximate its solution using a splitting method and computing 
characteristics. 

The first option is used in the numerical experiments. 

Remark 1. Numerical evidence reveals that using only one projection per stage 
(namely the one in the term haj jUB'^AYj) is sufficient to obtain up to second order 
of convergence. In this special case all the projections within the exponentials 
are removed, and the exponentials are simply treated as explained in the next 
section projDIRKCFNS:sec. 



3.4. Projected DIRK-CF methods for Navier-Stokes equations 

In this section we present an alternative semi-Lagrangian approach compared 



to the previous section 3.3 Also in this case the integration methods are a variant 



of the exponential integrators methods of |i7J. 
We rewrite ([19]) in the form 

y = B-^Ay - HB'^A + C(y))y + B-^C(y)y. (25) 

This is a differential equation on the subspace of discrete divergence-free vector 
fields, i.e. Dy = for all t. To approximate the solution of this equation, we 
consider a projection method of the type reviewed in [[241 IV. 4], see also [[T3l 
Sect. 5. 3. 3] and f25', Sect. VII. 2]. The idea is to use a one-step integrator 0/; for 



advancing the numerical solution of (|25j) by one step, and an orthogonal projection 
on the subspace of divergence-free vector fields applying 11 at the end of each step. 
We choose (ph to be the following integration method, in which the coefficients of 
both the DIRK-CF method and the underlying IMEX method are used: 

• the term (/ - H)B'^A y is treated implicitly with the DIRK coefficients. 
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• the term HB ^C{y)y is treated explicitly with the coefficients of the 
underlying explicit method, 

• the term B~^C{y)y is treated with the coefficients of the corresponding CF 
method. 

The projection 11 is used to guarantee divergence-free numerical approximations, 
i.e. = 0. We obtain 

for z = 1 : 5^ + 1 do 

Yi = iPiJn + h 2 ipiip-.' (ai,jIiB''AYj - aijHB-'C{Yj)Yj) + haijUB-'AYi 

7=1 

Yj:= j:kalyYkiovy=l,...,J 
cpi = exp(/75-iC(y/)) • ■ ■ QX^{hB-^C{Yl)) 
end for 

where 

f^s+ij = bj, a^+ij = bj, CK^+i^ = /3y, j,k = I, . . . ,s, = 0, a.s+i_i+i = 0. 

Since 11 is an orthogonal projection, the order of the method is not affected 
by the use of the projection in the update step|^ This approach requires only 
exponentials of pure convection problems, this will ease the implementation of 
the method as a semi-Lagrangian method. The exponentials 

exp(/i5~^C(w)) Wo 



^The projection does not need to be orthogonal, but should be guaranteed not to compromise 
the order of the method, an orthogonal projection will have this property. The target of the 
orthogonal projection map on the discrete divergence-free subspace is the element of shortest 
distance to the point which is projected. Since yn+i = TlYs and HyCfn+i) - Ys\\ = Oih''^^), where r is 
the order of the integration method, then 

IWf„+i) -y^+ill < llXfn+i) - + \\Y,-y„^,\\ = 0{h'-^\ 

because 

\\yn^A-ys\\<\\y{tn^l)-Y,\\. 
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correspond to pure convection problems: 



Uj = (w • V) u 

(26) 

u(0, x) = Wo(x) 

and have solutions 

u{t, X) = Wo(^(0), ^ = w(^), ^(0) = X. 

Here w(x) and Wo(x) are the interpolated functions obtained from w and wo, the 
intermediate numerical approximations of the velocity field. 

Observe that at each stage F, does not necessarily satisfy DF, = 0. 

4. Implementation issues 

Before proceeding to the numerical experiments, we describe some of 
the implementation issues, related to the discretization of the Navier-Stokes 
equations, and to the use of spectral element methods. 

4.1. Pressure-splitting scheme 

This scheme is used to obtain a cost efficient computation of solutions of 
discrete linear Stokes systems (see e.g.,[fT7]). 

The IMEX and DIRK-CF methods described in sections 3.2 and 3.3 give rise 
to linear Stokes systems of the form 

- D^Zi = Bfi 

(27) 

DYi =gi 

at each stage /, where 'H = ^5 - A, while fi, gi incorporate the vector fields at 
earlier stage values F/for j < i) and the contributions at boundary nodes (at stage 
i). 

For a method of order 1 or 2, the pressure-splitting scheme can be carried out 
in the following stepsj^ 

Stepl: 'HYi- D^zn = Bfi 

Step 2: DB-'D^Azi = -^.(DYi - g^) 

Step 3: Yi = Yi + Uijh B'^D^Azi, Zi = Zn + ^Zi- 



^'Extension to higher order methods is straightforward (See |29i and references therein). 
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The first step is an explicit approximation of the stage value of the velocity using 
the initial pressure z„. This approximation is not divergence-free. Steps 2 and 3 
are thus the projection steps which enforce the algebraic constrain and correct the 
velocity and pressure. Note that this approximation introduces a truncation error 
of order 0{h^), and is thus sufficient for methods of order up to 2 (see e.g.fTTll). 
Solving pT] ) directly would lead to solving linear equations with the operator 
]yH~^D^ for the Z,. However, the cost of inverting DI-f'^D^ is much higher than 
for inverting DB'^D^ in Step 2, since B is usually diagonal or tridiagonal and 
easier to invert than (which is usually less sparse). This explains the main 
advantage for using the pressure-splitting schemes in the numerical computations. 
We have exploited this advantage in the numerical experiments presented in 
sections 15.31 and 15. 4[ 

The matrix is a discrete Helmholtz operator and is symmetric positive- 
definite (SPD); the mass matrix B is diagonal and SPD, and thus easy to invert. We 
use conjugate gradient methods for "K ^ with B^ taken as preconditioner. The 
entire system pTj ) forms a symmetric saddle system, which has a unique solution 
for Yi provided D is of full rank. The choice of spatial discretization method 
guarantees this requirement. The system can be solved by a Schur-complement 
approach (or block LU-factorization) and the pressure-splitting scheme. 

Finally, we remark that the use of the pressure splitting scheme with our 
methods leads to overall approaches which can be regarded as a conventional 
projection schemes in the sense of [23] . 

4.2. Boundary conditions and discrete stiffness summation 

For the sake of completeness, we illustrate the strategy for implementing the 
boundary conditions in the context of spectral element methods. We use the 
spectral element notion known as the direct- stifl'ness summation (DSS), see for 
instance [|T2ll . 

Suppose we have to impose periodic or homogeneous Dirichlet boundary 
conditions, and that the variable y represents the values of the numerical solution 
at all discretization nodes in the computational domain (including boundary 
nodes). The variable y represents the restriction of y to the minimum degrees 
of freedom k required to define the numerical solution, while y contains typically 
redundant components. Thus if the number of components of y is N, then k < N . 
We denote by 2 a prolongation or "scatter" operator such that y = Qy. Associated 
to 2 is a restriction or "gather" operator denoted by . The operator Q is a 
N X k constant matrix of rank k < N . The variable y is referred to as the local 
variable, while y is the global variable. The DSS operator QQ^ ensures inter- 
element continuity and the fulfillment of the appropriate boundary conditions. 
So, for example, if the boundary conditions are periodic, given a vector y in the 
solution space or in the space of vector fields, QQ^y is periodic. 
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The spectral element discretization of the Navier-Stokes equations yields the 
discrete system 



B^=Ay + C{y)y - D^z (28) 
Dy = 0, (29) 

where y is assumed to be in the range of Q (i.e. y = Qy). The relation between the 
local and global operators is B = B Q, A = AQ, C(y) = C(Qy) Q and 
D = DQ. 

Applying on both sides of ([28]) by the DSS operator QQ^ we obtain 



= QQ^Ay + QQ^Cm - QQ^D^z, (30) 
Dy = 0. (31) 

The matrix S = QQ^B is N x N and invertible on the range of Q. In practice 
the integration methods are reformulated for the local variable y and the local 
operators. 

Indeed, in the computations the full data for the local variable y is stored, since 
all computations involving the operators B ^l-I'^ and {DB'^D^y^ must be done 
within the range of Q. These operators are symmetric and positive-definite, and 
so they can be inverted using a fast iterative solver such as the conjugate gradient 
method. For example the problem y = 'H~^ f is reformulated as follows: Find y 
such that 

\ Qy = y 

where 'ff = -^B - A. 

We refer to [flTll for further details on DSS and boundary conditions. In the 
experiments reported in this paper, no special treatment has been taken to enforce 
pressure boundary conditions, since the discrete pressure space is not explicitly 
defined on discretization nodes on the boundary. 

4.3. Reformulation of the integration methods in local variables 

In this section we briefly discuss the correct implementation of the methods of 
section 3|as applied to PU] ). The purpose of reporting here these implementation 



details is to explicitly highlight when care has to be taken in the implementation. 
See also the remark below. 

We first perform the elimination of the discrete pressure Lagrangian multiplier 



from ( [30| ) in analogy to ( [18] ). We use 

H = D^iDB-^DT^DB^^Q^, 
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and we get a system of ODEs for the variable y: 

^y = QQ^Ay + QQ^dy) y - QQ^H(A y + C(y) y), (33) 

Introducing the projection II = I-H allows us to write the following projected 
system of ODEs 

2y = QQ^flAy + QQ^flCm- (34) 



Denote 11 = 2115 Q . Applying the method of section 3.2 to the ordinary 



differential equation ([34]) split in its projected convection and projected diffusion 



terms, and then rewriting it as a method for the differential-algebraic equation ( [30| ) 
we obtain 

for ?■ = 1 : 5 do 

(S - haijQQ^A) Yi + hauQQ'D^Zi = Sy^ + /z 2 fl 2 UjA + CiYj)) Yj 
DYi = 
end for 

yn+i = Ys, 

under the assumption yn = Qyn- Since yn+i = Y^, the approximation of the velocity 
satisfies the discrete incompressibility constraint, Dy^+i = 0. 

Remark 2. We observe that this is not equivalent to what we obtain applying 



directly the IMEX method to (30), which written in ODE form is 



= QQ^Ay + QQ^Cm - QQ'HiAy + Cm)- 

In fact if we apply the IMEX method to ( [SO] ) we need to treat the term QQ^D^z 
either with the implicit method or with the explicit method, while in the approach 
outlined in this section we treated QQ^ HC{y)y explicitly and QQ^HAy implicitly, 
see also [I28H . 



Analogously, the method of section [33] applied to for ( 134] ) becomes 
for / = 1 : 5 + 1 do 

(2 - hauQQ^A)Yi + ha^QQ^ D'^Zi = i:Q^piy„ + h 2 a-.j!. Qip-.^p-.^ AYj 

7=1 

DYi = 

Yj:= i:,4/,for7=l,...,7 
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iPi = exp(m5-iC(y/)) • ■ ■ &x^{hIlB-^C{Y])) 
end for 

Jn+l = Ys+l- 



The method of section 3.4 may be reformulated in a similar way. 



5. Numerical experiments 

For the numerical experiments we shall employ a spectral element method 
(SEM) based on the standard Galerkin weak formulation as detailed out in [jTSl. 
We use a rectangular domain consisting of = Nx'x Ny uniform elements. The 



resulting discrete system has the form (17) or (30) (see section 3.1). The semi 



Lagrangian schemes associated to all the DIRK-CF methods in this section are 
achieved by tracing characteristics and interpolating as in [21]. The fourth order 
explicit RK method is used for approximating the paths of characteristic. 

5. 1. Temporal order tests for the IMEX methods 

We investigate numerically the temporal order of convergence of some 



IMEX methods following the algorithm described in section 3.2 The methods 
considered here are the second and third order IMEX-RK schemes with stiffly- 
accurate and L-stable DIRK parts [HI. We refer to them as IMEX2L and IMEX3L 



respectively. They are given by the Butcher tableaus in Table 1 and Table 2 



Table 1: IMEX2L: y = (2 - V2)/2 and 5 = 1 - 1 /(2y) 
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In the first example we consider the Taylor vortex problem with exact solution 
given by 

' u\ = - cos(7TXi) sm(nx2) e^p(-2n^t /Re), 
U2 = sm{nxi)cos(nx2)ei^p{-2n^t/Re), 



1 2 
--[cos(2;r.xi) + cos(2;r;v:2)] exp(-47r t/Re), 



(35) 



where Re = 1/v is the Reynolds number, and u := {ui,U2), x := (xi,X2). The 
boundary condition is doubly-periodic on the domain Xi, X2 6 [-1,1], and we 
choose Re = In^. The initial conditions are determined from the exact solution 
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Table 2: IMEX3L 
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p5| ). For the spatial discretization we use a spectral method of high order N = \2, 
with Ne = 4 elements, and the time integration is done up to time T = I. For each 
time-step h = T k = \, . . . ,9, the global error between the numerical solution 
and the exact PDF solution (at time T) are measured in the Hi- and L2-norm^ 
respectively, for the velocity and pressure. These are illustrated in log-log plots of 
the errors against the time-steps. The results for both the IMFX2L and IMFX3L 



show temporal convergence of order 2 and 3 respectively (see Figure 1 



5.2. Temporal order tests for the DIRK-CF methods 

Using the IMFX2L and IMFX3L methods, we construct two DIRK-CF 
methods, namely, DIRK-CF2L and DIRK-CF3L, of classical orders 2 and 3 



respectively. Both DIRK-CF methods are applied to ( 17) following the algorithm 



discussed in section 3.3[ To approximate the exponentials, we use a semi- 



Lagrangian approach coupled with a high order approximation based on ([24]) with 
E X E2OT £3 where 

/j2 /j3 

E2 = -—B-^CHB-^C and E^ = E2 + —B-^C[(nB-^Cf - (B-^Cf] 
2 6 

for order 2 and order 3 DIRK-CF methods respectively. These are error terms of 
order 0{h^) and 0{h^) obtained via Taylor series expansions of the exponentials 
in ( |24l ). The exponential on the right hand side of (|24]) is accurately computed 
using semi-Lagrangian methods for pure convection problems of the type (|26]). 
Using the same test example as for section |5]2 we observe the temporal order of 



convergence 2 and 3, in both the velocity and pressure (see Figure 1| ). 



Following the semi-Lagrangian algorithm of section 3.4 we apply the second 



*See 



Appendix A. 2 



; for the definitions of these norms. 
^ The slope is measured by a linear fitting over portions where the temporal errors are dominant 
over the spatial errors, in the given log-log plot. The reported value of the slope indicates the 
approximate temporal order of convergence. 
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Figure 1: Temporal order of convergence. Test problem: Taylor vortex (35 1; Re = 2n^, 
T ^ \, N = 12, Ne = 4, = Ny = 2, Q = [-1,1]^, h = Tjl^, k = 1,...,9. be: 
periodic, (a) velocity error (Hi): IMEX2L (slope = 2.0084), IMEX3L (slope = 2.9912); 
(b) pressure error (L2): IMEX2L (slope = 1.9097), IMEX3L (slope - 2.9035); (c) 
velocity error (Hi): DIRK-CF2L (slope = 2.1485), DIRK-CF3L (slope = 2.9212); (d) 
pressure error (L2): DIRK-CF2L (slope - 2.0366), DIRK-CF3L (slope - 2.7313). 



and third order DIRK-CF methods, this time for the test problem [[22| with exact 
solution given by 

Ml = nsm(2nx2) sm^(nxi) sm(t), 
' U2 = -7Tsm{2n xi)sm^(7T X2)sm(t), (36) 
p = cos(;r.x:i) sin(;rjC2) sin(0, 

for Xi, X2 6 [0, 1] and t e [0, T], with T = I. A corresponding forcing term f is 
added to the momentum equations ([T]) so that ( [361 ) the exact solution. In this test 
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Figure 2: Temporal order of convergence. Test problem (36i; Re - 100, T - I, N 



8, Ne = 16, A^., 



4,a = [0, 1]^ h = h = TIT, it ^ 5, ... ,9. be: homogeneous 



Dirichlet. (a) velocity error (L2): DIRK-CF2L (slope - 2.1268), DIRK-CF3L (slope 
= 2.8632); (b) pressure error (L2): DIRK-CF2L (slope - 2.1296), DIRK-CF3L (slope 
= 3.1628). 



case we have used Re = 100. Meanwhile (|36]) is used to prescribe the initial data 
and boundary conditions (homogeneous Dirichlet on the entire boundary). The 
errors are all measured in the L2-norm. We observe temporal order of convergence 
2 and 3, in both the velocity and pressure (see Figure 2[ ). 

In the subsequent sections [53 and [54 we present a set of numerical 
experiments that illustrate the potentials of the semi-Lagragian exponential 
integrators [7] for the treatment of convection-dominated problems. Two 
examples involving the incompressible Navier-Stokes models at high Reynolds 
numbers are considered. These examples are the shear-layer roll up problem 
in [in [161 |T8]|, and the 2D lid-driven cavity problem (see Il20l [4ll and references 
therein). The second order semi-Lagrangian DIRK-CF2L method (named SL2L 
in [[3) is used in each of these experiments. The pressure-splitting technique [[TTII 



(discussed in section 4.1) is applied to solve the discrete linear Stokes system 
that arises at each stage of the DIRK-CF method. The results reported in both 
sections [5.3[ and [5.4[ indicate that the semi-Lagrangian exponential integrators 



permit the use of large time-steps and Courant numbers. 

5.3. Lid-driven cavity flow in 2D 

We consider the 2D lid-driven cavity problem on a domain (jc, j) e Q := [0, 1]^ 
with initial data u = (w, v) = (0, 0) and constant Dirichlet boundary conditions 



1 1 on upper portion of dQ. 
1 elsewhere on dQ. 



V = on dQ.. 



(37) 
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We demonstrate the performance of the second order DIRK-CF method 
(SL2L, by the nomenclature of [|7]|). Spectral element method on a unit square 
domain [0, 1]^ with A^^ = 10 x 10 uniform rectangular elements and polynomial 
degree = 10 is used (see [36J). A constant time-step, h = 0.03, is used, 
corresponding to a Courant number of Cr « 9.0911. The time integration is 



carried out until the solution attains steady-state. The results in Figure 3 show 
the evolution of the center velocity (at Re = 400) up to steady state. It can be 
observed from this figure that steady state is attained at time t « 40. At steady state 
the relative error (L2-norm), between the velocity at a given time (?„+i) relative to 
the velocity at the preceding time (?„), has decreased to 0(lO'^). The results also 
match with those of [|36l. In Figure 4i-b we plot the streamline contours of the 
stream functions, choosing contour levels as in flU. Meanwhile in Figure 4:-d 
plots of the centerline velocities (continuous line, for Re = 400, dashed line, for 
Re = 3200) show a good match with those reported in [i20il (plotted in red circles). 



5.4. Shear-layer roll up problem 

We now consider the shear-layer problem 
with initial data u = (u, v) given by 



on a domain Q := [0, 1]^ 



u = 



ftanh(p(3;-0.25)) 
ltanh(p(0.75 -y)) 



for y < 0.5 
for y > 0.5 



V = 0.05 sin(2;r x) 



(38) 



which corresponds to a layer of thickness (9(1 /p). Doubly-periodic boundary 
conditions are applied. 





(a) u at center 



(b) V at center 



Figure 3: Results of the second order semi-Lagrangian DIRK-CF method (SL2L) for the 
2D lid-driven cavity problem. We have {x,y) e [0, 1]^; A^^ ^ 10 x 10, N = 10, h = 0.03, 
Cr = 9.0911, /?e = 400. (a) Evolution of the horizontal velocity component u at the 
domain center (x = 0.5, y = 0.5): t e (0, 112.08), (b) Evolution of the vertical velocity 
component v at the domain center (x = 0.5, y = 0.5): t € (0, 112.08). 
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(a) stream function {Re - 400) (b) stream function {Re = 3200) 




(c) center line velocity u (d) center line velocity v 



Figure 4: Results of a second order DIRK-CF method for the 2D lid-driven cavity 
problem. We have {x,y) e [0, 1]^; A^^ = 10 x 10, N = 10, h ^ 0.03, Cr = 9.0911. In 
blue continuous line (our numerical solution); in red circles (O, reference solution from 
EOll ). (a) Streamline contours of the solution for Re = 400, (b) Streamline contours of the 
solution for Re = 3200, (c) Horizontal velocity component u along the vertical center line 
{x = 0.5), (d) Vertical velocity component v along the horizontal center line (y = 0.5). 



In Figure 5 we demonstrate the performance of various second order methods 
including two DIRK-CF methods (SL2 & SL2L, by the nomenclature of [7l), 
and also a second order semi-Lagrangian multistep exponential integrator (named 
BDF2-CF2, in S). The results are obtained at time t = 1.5, using a filter-based 
spectral element method (see IfTSlI ) with A^g = 16 x 16 elements and polynomial 
degree N = S. The specified Reynolds number is Re = 10^, while p = 30 and 
time-steps used are h = 0.002, 0.005, 0.01 corresponding to a Courant numbers 
of Cr « 0.6393, 1.5981, 3.1963 respectively. The filtering parameter used in 
each experiment is or = 0.3 (see for example IITSlI). However, the time-step and 
Courant number are up to about 10 times larger than that report in [fTSl . The initial 
values for the BDF2-CF are computed accurately using the second order DIRK- 
CF (SL2L) with smaller steps. The results are qualitatively comparable with those 
in[[l6l[T8l|. 

In [Figure 6] we demonstrate the performance of the second order DIRK- 
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(a) BDF2CF: h = 0.002 (b) SL2L: h - 0.002 (c) SL2: h - 0.002 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



(d) BDF2CF: h - 0.005 (e) BDF2CF: h = 0.01 (f) SL2L: h - 0.01 

Figure 5: Results of second order DIRK-CF methods (SL2 & SL2L) and BDF2-CF 
method for the shear-layer roUup problem. We have ix,y) € [0, 1]^; Ne = 16 x 16 = 
256, N = S. (filtering, a = 0.3), p = 30, Re = 10^. Vorticity contours (-70 to 70 by 15) of 
the solution at time t = 1.5. The corresponding Courant numbers are (a) Cr = 0.6393, (b) 
Cr - 0.6393, (c) Cr - 0.6393, (d) Cr - 1.5981, (e) Cr = 3.1963, (f) Cr - 3.1963. 



CF method (SL2L). The results are obtained at times t = 0.8, 1.0, 1.2 and 1.5 
respectively, using spectral element method (without filtering) with A^g = 16 x 16 
elements and polynomial degree = 16. The specified Reynolds number is 
Re = 10^, while p = 30. The time-step used is h = 0.01, corresponding to a 
Courant number of Cr « 11.9250. This time-step is 10 times larger than that 
reported in [ITSll . Again the results are well comparable to those in [16, 18 1. 



Finally in Figure 7 we demonstrate the performance of the second order 
DIRK-CF method (SL2L) for the "thin" shear-layer roll up problem, so defined 
for p = 100. The results are obtained at times t = 0.8, 1 .0, 1 .2 and 1 .5 respectively, 
using spectral element method (without filtering) with A^g = 16 x 16 elements and 
polynomial degree = 16. The specified Reynolds number is i?e = 4 x 10"^. The 
time-step used is h = 0.01, corresponding to a Courant number of Cr ^ 1 1.9250. 
The results are well comparable to those in [fT6l[T^ , except that we used 10 times 
the step size in time. 
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(a)f-0.8 (b)f=1.0 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

X X 



{c)t=\.2 (d)f=1.5 

Figure 6: Results of second order DIRK-CF method (SL2L) for the shear-layer roUup 
problem. We have {x,y) e [0,1]^; A^^ - 16 x 16 = 256, N = \6, h ^ 0.01, Cr = 
1 1.9250, p = 30, Re = 10^. Vorticity contours (-70 to 70 by 15) of the solution at time (a) 
t = 0.8, (b) t = 1.0, (c) t - 1.2, (d) t - 1.5. 



6. Conclusion 

In this paper we have presented a class of semi-Lagrangian methods for the 
incompressible Navier-Stokes equations of high order in time to be used with 
high order space discretizations, such as for example spectral element methods. 
We have proposed a strategy to maintain the high temporal order also in the 
presence of constraints. As a by product, we have also derived projection methods 
based on IMEX Runge-Kutta schemes which have been used for comparison. The 
methods have been implemented and tested, and have been shown the predicted 
order of convergence in the case of periodic and no- slip boundary conditions. 
For convection-dominated test problems, in 2D with high Reynolds number, 
the semi-Lagrangian methods showed improved performance compared to their 
Eulerian counterparts, allowing for the use of considerably larger time-steps. So 
far methods up to order three have been implemented. Order four methods where 
obtained in [[8]|, and will be implemented and tested in future work. 
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(a) f - (b) t = 0.8 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

X X 



(c)f=1.0 (d)f=1.2 

Figure 7: Results of second order DIRK-CF method (SL2L) for the "thin" shear-layer 
rollup problem. We have {x,y) € [0, 1]^; A^^ ^ 16 x 16 = 256, N ^ 16, h ^ 0.01, Cr = 
1 1.9250. (no filtering), p = 100, Re = 40, 000. Vorticity contours (-36 to 36 by 13) of the 
solution at time (a) t = 0,{h)t = 0.8, (c) t = 1.0, (d) t = 1.2. 



Appendix 

Appendix A.l. IMEX and DIRK-CF methods for convection-diffusion equations 

A DIRK-CF method is a IMEX Runge-Kutta method of exponential type. It 
is defined by two Runge-Kutta tableaus one of implicit type to treat the linear 
diffusion and one of explicit type to treat the nonlinear convection by composition 
of exponentials. These methods have a Runge-Kutta like format with two sets of 
parameters: 

= {aij}ij=i,...,s, h = [bi,...,b,], c = [ci,...,cj 

and 

ajj, fi'i, i = I,. . . s, j = I,. . .,s, I = \,. . . ,J, c = [ci,...,cj. 
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When applied to the convection-diffusion problem 

U{t) + C{U{t)) U{t) = AU{t), 

with linear diffusion and nonlinear convection the DIRK-CF methods have the 
following format: 

for i=l:s do 

U = 'PiUn + hiipij(aijA^j) 

7=1 

ifi = cxpih 2^ aljCm)) ■ ■ ■ exp(/z 2, alcm)) 
(fij := (fjifij^ 
end for 

cp,^, = exp(/i Y.kP)C{nAu)) ■ ■ ■ exp(/z Y.kP\Cm)) 
The methods are associated to the two Butcher tableaus. 



(A.l) 



c 




c 






b ' 


A 

b 



where we have defined 



1=1 



1=1 



(A.2) 



for z = I,. . .,s, Ji = {a,,;},j=i,...,i and b = [bi, . . . The coefficients of the 
first tableau are used for the linear vector field Ay while the coefficients of the 



second tableau, split up in the sums ( |A.2[ ), are used for the nonlinear vector field 
C(y)y. We choose the first tableau to be a DIRK (diagonally implicit Runge-Kutta) 
method, this means we are solving only one linear system per stage. The tableaus 



( |A.1[ ) are typically chosen so that they define a classical IMEX method, which we 
call the underlying IMEX method, see [IJ and [7| for more details. This IMEX 
method has the format 



for / = 1 : 5 do 
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end for 

Un-,1 = Un + h Zt=i(biCmm + biAU). 



The order theory for classical IMEX methods reduces to the theory of 
partitioned Runge-Kutta methods, [24J. Given an implicit and an explicit method 
of order k they must satisfy extra compatibility conditions in order for the 
corresponding IMEX method to have order k. The extension of this theory to 
the DIRK-CF methods has been discussed in [|71 and [[8]|. 

Appendix A.2. Definition of norms 

For a square-integrable (respectively Hi) function u : Q ^ R", where Q. c W 
is bounded and connected, the L2-norm (|| • \\l2{si)) and the //i-norm (|| • Wniisi)) are 
defined by 



m\L2{ii) ■-- 



V (=1 
/ n 



1/2 



(uj + Vui ■ Vui) dQ 



1/2 



(A.3) 
(A.4) 



In the spectral element approximations the continuous integrals of numerical 
solutions are accurately computed using Gauss quadrature rules. 
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